function [x,zo]=estnoiseg(yf,tz,pp)
%ESTNOISEG - estimate MMSE noise spectrum [x,zo]=(yf,tz,pp)
%
% Usage:    ninc=round(0.016*fs);   % frame increment [fs=sample frequency]
%           ovf=2;                  % overlap factor
%           f=rfft(enframe(s,hanning(ovf*ninc,'periodic'),ninc),ovf*ninc,2);
%           f=f.*conj(f);           % convert to power spectrum
%           x=estnoiseg(f,ninc/fs); % estimate the noise power spectrum
%
% Inputs:
%   yf      input power spectra (one row per frame)
%   tz      frame increment in seconds
%           Alternatively, the input state from a previous call (see below)
%   pp      algorithm parameters [optional]
%
% Outputs:
%   x       estimated noise power spectra (one row per frame)
%   zo      output state
%
% The algorithm parameters are defined in reference [1] from which equation
% numbers are given in parentheses. They are as follows:
%
%        pp.tax      % smoothing time constant for noise power estimate [0.0717 seconds](8)
%        pp.tap      % smoothing time constant for smoothed speech prob [0.152 seconds](23)
%        pp.psthr    % threshold for smoothed speech probability [0.99] (24)
%        pp.pnsaf    % noise probability safety value [0.01] (24)
%        pp.pspri    % prior speech probability [0.5] (18)
%        pp.asnr     % active SNR in dB [15] (18)
%        pp.psini    % initial speech probability [0.5] (23)
%        pp.tavini   % assumed speech absent time at start [0.064 seconds]
%
% If convenient, you can call estnoiseg in chunks of arbitrary size. Thus the following are equivalent:
%
%                   (a) dp=estnoiseg(yp(1:300),tinc);
%
%                   (b) [dp(1:100,:),z]=estnoiseg(yp(1:100,:),tinc);
%                       [dp(101:200,:),z]=estnoiseg(yp(101:200,:),z);
%                       [dp(201:300,:),z]=estnoiseg(yp(201:300,:),z);


% This is intended to be a precise implementation of [1] for a frame rate of 62.5 Hz.
% Time constants are adjusted for other frame rates.
%
% Refs:
%    [1] Gerkmann, T. & Hendriks, R. C.
%        Unbiased MMSE-Based Noise Power Estimation With Low Complexity and Low Tracking Delay
%        IEEE Trans Audio, Speech, Language Processing, 2012, 20, 1383-1393

%	   Copyright (C) Mike Brookes 2012
%      Version: $Id: estnoiseg.m 2268 2012-09-09 11:40:59Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nr,nrf]=size(yf);          % number of frames and freq bins
x=zeros(nr,nrf);            % initialize output arrays
if isempty(yf) && isstruct(tz)             % no real data
    zo=tz;                  % just keep the same state
else
    if isstruct(tz)       % take parameters from a previous call
        nrcum=tz.nrcum;   % cumulative number of frames
        xt=tz.xt;          % smoothed power spectrum
        pslp=tz.pslp;               % correction factor (9)
        tinc=tz.tinc;     % frame increment
        qq=tz.qq;         % parameter structure
    else
        tinc = tz;          % second argument is frame increment
        nrcum=0;            % no frames so far
        % default algorithm constants
        qq.tax=0.0717;      % noise output smoothing time constant = -tinc/log(0.8) (8)
        qq.tap=0.152;       % speech prob smoothing time constant = -tinc/log(0.9) (23)
        qq.psthr=0.99;      % threshold for smoothed speech probability [0.99] (24)
        qq.pnsaf=0.01;      % noise probability safety value [0.01] (24)
        qq.pspri=0.5;       % prior speech probability [0.5] (18)
        qq.asnr=15;         % active SNR in dB [15] (18)
        qq.psini=0.5;       % initial speech probability [0.5] (23)
        qq.tavini=0.064;        % assumed speech absent time at start [64 ms]

        if nargin>=3 && ~isempty(pp)  % update fields from pp input
            qqn=fieldnames(qq);
            for i=1:length(qqn)
                if isfield(pp,qqn{i})
                    qq.(qqn{i})=pp.(qqn{i});
                end
            end
        end
    end

    % unpack parameter structure

    tax=qq.tax;    % noise output smoothing time constant = -tinc/log(0.8) (8)
    tap=qq.tap;    % noise output smoothing time constant = -tinc/log(0.9) (23)
    psthr=qq.psthr;  % threshold for smoothed speech probability [0.99] (24)
    pnsaf=qq.pnsaf;      % noise probability safety value [0.01] (24)
    pspri=qq.pspri;       % prior speech probability [0.5] (18)
    asnr=qq.asnr;         % active SNR in dB [15] (18)
    psini=qq.psini;       % initial speech probability [0.5] (23)
    tavini=qq.tavini;     % assumed speech absent time at start [0 seconds]

    % derived algorithm constants

    ax=exp(-tinc/tax); % noise output smoothing factor = 0.8 (8)
    axc=1-ax;
    ap=exp(-tinc/tap); % noise output smoothing factor = 0.9 (23)
    apc=1-ap;
    xih1=10^(asnr/10); % speech-present SNR
    xih1r=1/(1+xih1)-1;
    pfac=(1/pspri-1)*(1+xih1); % p(noise)/p(speech) (18)

    if ~nrcum       % initialize values for first frame
        xt=psini*mean(yf(1:max(1,min(nr,round(1+tavini/tinc))),:),1);  % initial noise estimate
        pslp=repmat(psini,1,nrf); % initialize smoothed speech presence prob
    end

    % loop for each frame
    for t=1:nr
        yft=yf(t,:);        % noisy speech power spectrum
        ph1y=(1+pfac*exp(xih1r*yft./xt)).^(-1); % a-posteriori speech presence prob (18)
        pslp=ap*pslp+apc*ph1y; % smoothed speech presence prob (23)
        ph1y=min(ph1y,1-pnsaf*(pslp>psthr)); % limit ph1y (24)
        xtr=(1-ph1y).*yft+ph1y.*xt; % estimated raw noise spectrum (22)
        xt=ax*xt+axc*xtr;  % smooth the noise estimate (8)
        x(t,:)=xt;  % save the noise estimate
    end
end
if nargout>1    % we need to store the state for next time
    zo.nrcum=nrcum+nr;      % number of frames so far
    zo.xt=xt;          % smoothed power spectrum
    zo.pslp=pslp;               % correction factor (9)
    zo.tinc=tinc;     % must be the last one
    zo.qq=qq;
end
if ~nargout
    clf;
    subplot(212);
    plot((1:nr)*tinc,10*log10([sum(yf,2) sum(x,2)]))
    ylabel('Frame Energy (dB)');
    xlabel(sprintf('Time (s)   [%d ms frame incr]',round(tinc*1000)));
    axisenlarge([-1 -1.05]);
    legend('input','noise','Location','Best');
    subplot(211);
    plot(1:nrf,10*log10([sum(yf,1)'/nr sum(x,1)'/nr]))
    ylabel('Power (dB)');
    xlabel('Frequency bin');
    axisenlarge([-1 -1.05]);
    legend('input','noise','Location','Best');
end

